Direct calculation of the hard-sphere crystal/melt interfacial free energy 
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We present a direct calculation by molecular-dynamics computer simulation of the crystal/melt 
interfacial free energy, 7, for a system of hard spheres of diameter a. The calculation is performed by 
thermodynamic integration along a reversible path defined by cleaving, using specially constructed 
movable hard-sphere walls, separate bulk crystal and fluid systems, which are then merged to form an 
interface. We find the interfacial free energy to be slightly anisotropic with 7 = 0.62±0.01, 0.64±0.01 
and 0.58±0.01fcflT/(T^ for the (100), (110) and (111) fee crystal/fluid interfaces, respectively. These 
values are consistent with earlier density functional calculations and recent experiments measuring 
the crystal nucleation rates from colloidal fluids of polystyrene spheres that have been interpreted 
[Marr and Gast, Langmuir 10, 1348 (1994)] to give an estimate of 7 for the hard-sphere system of 
0.55 ± 0.02fcsr/cr^, slightly lower than the directly determined value reported here. 

PACS number(s): 68.45-v, 05.70.Np, 05.10.-a, 68.35.Md 
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A detailed microscopic description of the interface be- 
tween a crystal and its melt is necessary for a full under- 
standing of such important phenomena as homogeneous 
nucleation and crystal growth Computer simula- 

tion studies of model materials have had some success 
in elucidating the phenomenology of such systems [Q, 
the importance of such work being enhanced by the near 
absence of reliable experimental studies. These efforts, 
however, have been primarily focused on structural and 
dynamical properties, since the central thermodynamic 
property, the crystal/melt interfacial free energy, is dif- 
ficult to measure by simulation or experiment. In this 
work we report the results of a direct calculation via 
molecular-dynamics (MD) simulation of the crystal/melt 
surface free energy of the hard-sphere system, one of the 
most important reference models for simple materials. 

The crystal/melt surface free energy, 7, is defined as 
the (reversible) work required to form a unit area of in- 
terface between a crystal and its coexisting melt. Exper- 
imentally, 7 can be measured either indirectly from mea- 
surements of crystal nucleation rates interpreted through 
classical nucleation theory, or directly by contact angle 
measurements [0. Using the former method, TurnbuU 
1^ estimated 7 for a number of materials and found a 
strong empirical correlation between the values obtained 
and the latent heat of fusion for each material given by 
the relation 7 w Cr^fHp^^^, where p is the number den- 
sity of the crystal and with Ct (the TurnbuU coefficient) 
taking on the value 0.45 for most metals and 0.32 for 
other mostly nonmetallic materials. For the hard-sphere 
system considered in this work, recent experiments ^] 
of the crystallization kinetics of a colloidal suspension of 
uncharged polystyrene spheres, which closely mimic hard 
spheres, have been interpreted within a classical nucle- 
ation model to yield an estimate for 7 of the hard-sphere 
system of 0.55±0.02kBT/a^ 0. This value is in agree- 
ment with that predicted using the empirical rela- 
tionship above assuming a Ct of 0.45 and values of AfH 
and coexistence densities for hard spheres as determined 
by MD simulation M. Unfortunately, the accuracy of 



the values of 7 obtained from nucleation rates is severely 
limited by the approximations inherent in classical nu- 
cleation theory. More accurate values can be obtained 
directly using contact angles, but such experiments are 
difhcult and only a few materials have been studied to 
date. One notable example is a series of grain boundary 
contact angle experiments on bismuth that deter- 
mined 7 to be relatively independent of crystal orienta- 
tion at 61.3x10"'^ J/m^, which is about 10% higher than 
TurnbuU's estimate from nucleation rates of 54.4x10^'^ 
J/m2. 

In recent years, the primary theoretical approach to 
studying the structure and thermodynamics of the crys- 
tal/melt interface has been density- functional theory 
(DFT) pd|-p^. For these studies, the hard-sphere sys- 
tem has been the benchmark calculation, due to the sim- 
plicity of the interaction and the availability of accurate, 
analytical formulas for the properties of the fluid. How- 
ever, as discussed by Marr ||] the value of 7 obtained 
is highly dependent on the approximations used to con- 
struct the DFT and the reported values range from 0.25 
to i.OOksT/a^ . The DFT studies also disagree dramat- 
ically in the degree of orientation dependence of the in- 
terfacial free energy. Unfortunately, in the absence of 
simulation results it has been difficult to assess the valid- 
ity of the individual approaches, although only the DFT 
approach of Curtin ||l^ (7 = 0.62kBT/a'^) and the re- 
lated one of Marr and Gast |Q (7 = O.GOfcsT/cr^) come 
close to the nucleation result of 0.55 ± 0.02kBT/a-^ . 

To date, the only reliable calculation of the crys- 
tal/melt interfacial free energy via simulation is that of a 
system of particles interacting via a truncated Lennard- 
Jones potential by Broughton and Gilmer ||l^. In that 
work, a series of continuous, external "cleaving" poten- 
tials are used to separate (cleave) separate samples of 
bulk liquid and fee crystal, prepared at the calculated 
coexistence temperature and densities. The solid and 
liquid slabs thus produced are then placed next to one 
another and the cleaving potentials removed to merge 
them into a coexisting interface. The reversible work to 
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perform these steps can be calculated by thermodynamic 
integration, giving a direct calculation of 7 for this sys- 
tem. The values of 7 at the triple point were found to be 
statistically isotropic with ^(p- jt = 0.35±0.02, 0.34±0.02 
and 0.36 ±0.02 for the (111), (100) and (110) interfaces, 
respectively. 

For the hard-sphere system, the Broughton-Gilmcr 
procedure must be modified since the algorithm for MD 
simulation of discontinuous hard-core potentials is con- 
ceptually very different from the algorithm for continu- 
ous potentials. The latter are performed by integrating 
the system of ordinary differential equations, while in the 
former the dynamical algorithm proceeds on a collision 
by collision basis. Therefore, incorporating continuous 
cleaving potentials into the coUisional algorithm would 
result in lost efficiency and substantial modification of 
the structure of the algorithm. 

In this Letter we introduce an approach, which uses 
only hard-sphere interactions in order to cleave the 
bulk hard-sphere systems. This allows us to apply the 
Broughton-Gilmer cleaving procedure to the hard-sphere 
system with only minor changes to the algorithm struc- 
ture. The idea of our approach is illustrated in Fig. |l]. To 
cleave the bulk system at a cleaving 'plane (shown by the 
dashed line in Fig. |]), the spheres are assigned types 1 
or 2 based on their position relative to the plane. Next, 
two walls, which are also assigned types 1 and 2, are 
placed on the opposite sides of the cleaving plane. The 
two types are introduced in order to specify interaction 
between the spheres and the walls; namely, the walls in- 
teract only with the spheres of similar type. Therefore, 
when the walls are placed as shown in Fig. ^ and the 
distance from the walls to the cleaving plane is larger 
than the sphere radius, the walls do not interact with 
the spheres. It is important that, during a simulation 
run, a sphere changes its type whenever it crosses the 
cleaving plane. Because of the periodic boundary condi- 
tions in the z direction, another plane must be defined 
sufficiently far away from the cleaving plane, where the 
spheres also change type. 

The cleaving of the system is achieved by slowly mov- 
ing the walls towards each other (as shown by the arrows 
in Fig. starting from the initial position Zi, where 
the walls do not interact with the system, till Zf, where 
the spheres of different types no longer collide with each 
other at the cleaving plane. During the process, the aver- 
age pressure on the walls, P{z), is measured as a function 
of wall position. The work per unit area required to per- 
form the cleaving is given by the integral 

w = ^ ^ P{z) dz . (1) 

Thus the crystal-fluid interfacial free energy, 7, can be 
measured in the reversible process involving the follow- 
ing four steps: (1) Cleave the bulk crystal by inserting 
two walls at the cleaving plane and moving them from 
Zi to Zf] (2) Cleave the bulk fluid in a similar way; (3) 
Juxtapose the cleaved crystal and fluid systems by chang- 
ing the periodic boundary conditions while retaining the 



FIG. 1. Diagram illustrating the cleaving of the bulk 
hard-sphere system by two moving walls. Spheres are as- 
signed types 1 and 2 based on their position with respect to 
the cleaving plane (dashed line). Two walls of types 1 and 2, 
which interact only with spheres of similar type, are placed 
on the opposite sides of the cleaving plane, so that initially 
there are no collisions between walls and spheres (as shown 
on the diagram). The system is then cleaved by moving the 
walls in directions indicated by the arrows. 

crystal and fluid systems restricted by the respective 
cleaving walls; (4) Slowly move the walls back to their 
initial positions with respect to the cleaving planes. This 
series of steps is the same as that used by Broughton and 
Gilmer |p7| , except that, in our case, no work is done on 
the system in step 3. The interfacial free energy is given 
by the sum 7 = wi + W2 + u;4 , where W4 is negative and 
consists of the work done by the walls on the crystal and 
fluid parts of the system during the process of removing 
the cleaving walls. 

The structure of the cleaving walls is crucial to the 
success of the procedure. The main requirement is that 
the insertion of the walls perturbs the systems as little 
as possible. Our approach is to use walls made of layers 
of ideal crystal oriented in correspondence with the ori- 
entation of the crystal system. For the (100) and (111) 
orientations it is sufficient to use one layer, while for the 
(110) orientation we use two layers. Such a choice of 
the wall structure ensures minimal perturbation of the 
cleaved crystal, while the cleaved fiuid is expected to 
form properly oriented interfacial layers. On a techni- 
cal side, the implementation of such a wall structure is 
quite simple, since we can treat collisions with the walls 
in the same manner we treat collisions between all the 
spheres in the system, except that the spheres forming 
the walls are assigned infinite mass. 

The position of the walls, z, is measured by the dis- 
tance of the centers of the spheres forming the walls to 
the cleaving plane. Obviously, the walls do not interact 
with the system when z > a. Therefore, we set the ini- 
tial position of the walls at zi = a. The pressure P{z) 
is obtained by moving the walls from Zi to Zf in steps 
of O.Olcr. In order to move the walls to a new position, 
we assign a small velocity (about 0.1% of the average 
particle velocity) to the spheres forming the walls, and 
run the simulation until the walls reach the new posi- 
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FIG. 2. Step 1: Cleaving crystal system. Pressure on the 
cleaving walls as a function of wall position for the three ori- 
entations of the interface. The error bars are smaller than the 
size of the symbols. 



FIG. 3. Step 2: Cleaving fluid system. Pressure on the 
cleaving walls as a function of wall position for the three ori- 
entations of the interface. The error bars are of the order of 
the size of the symbols. 



tion. At that moment the wall velocity is set to zero, 
and the velocities of the particles are rescaled in order 
to restore the initial value of the total kinetic energy of 
the system. Before measuring the pressure, we allow the 
system to relax in an equilibration run. In order to ver- 
ify the reversibility of the cleaving process, we have also 
simulated the reverse process and measured the pressure 
while moving the walls from Zf back to Zi. The details 
of the simulation process and obtained results follow. 

Step 1: Cleaving the crystal For each of the three 
orientations, we start with a crystal at a density pc = 
1.037(T~^, which corresponds to the crystal-fluid coexis- 
tence pressure of 11. SSfesTcr"'^ [|l8|. In order to minimize 
the size effects and the amount of stress in the crystal 
introduced by the cleaving walls, we use large systems 
of about 8000 spheres and approximate dimensions of 
14tT X 14cr X 40(T. (To perform the simulations efhciently 
for such large systems, we use the algorithm of Rappaport 
.) The cleaving plane is placed in the middle between 
two crystal layers. The dependence of the pressure on the 
wall position is shown in Fig. ||. The walls do not inter- 
act with the crystal until they move sufficiently close to 
the crystal layers (about O.Tcr for all orientations). Then 
the pressure on the walls quickly rises to slightly above 
the bulk crystal pressure. The steepness of the rise is 
directly related to the compactness of the layers for each 
orientation. The final positions, Zf, where the spheres of 
different types no longer collide across the cleaving plane, 
were determined to be O.Slcr, 0.16cr, and 0.35cr for the 
(100), (110), and (111) system orientations, respectively. 
No hysteresis was observed in the reverse process. 

Step 2: Cleaving the fluid The fluid systems consisting 
of about 7400 particles are prepared at the coexistence 
density p/ = 0.939cr~^ using box dimensions nearly iden- 
tical to the crystal systems. Unlike in step 1, the cleaving 
walls begin to interact with the fluid system as long as 
z < a. As can be seen in Fig. H, the pressure on the walls 



increases approximately linearly until the fluid near the 
cleaving walls begins to develop signiflcant crystal-like 
ordering commensurate with the wall structure. At that 
point the pressure in the bulk fluid decreases to about 
11.2kBT(T~^ , after which the dependence of pressure on 
the wall position follows essentially the same curve as 
during the cleaving of the crystal system, which leads to 
the same values of as in step 1. 

Simulation of the reverse process shows that the or- 
dering of the fluid against the walls is the source of some 
hysteresis. However, we have found that the magnitude 
of the hysteresis can be always reduced to within the sta- 
tistical error by increasing the duration of the equilibra- 
tion run. In other words, at every position of the cleaving 
walls, the pressure on the walls eventually converges to 
the same value (within the statistical error bounds) in 
both forward and reverse processes. 

Step 3: Changing boundary conditions The combined 
system has two cleaving planes and four walls. The crys- 
tal part of the system and the two walls restraining it are 
assigned type 1, while the fluid part and the other two 
walls are assigned type 2. 

Step 4- Removing the cleaving walls As can be seen 
in Fig. ^, the pressure on the walls restraining crystal 
and fluid parts of the system is essentially the same as 
in steps 1 and 2, respectively, except that the fluid part 
retains its structure in the interfacial region. Thus the 
main contribution to the interfacial free energy comes 
from the pressure of the fluid on the cleaving walls before 
signiflcant crystal-like ordering at the wall develops. 

The work done during each step and the resulting in- 
terfacial free energy for each of the three orientations is 
given in Table ||. The average values of 7 for the three ori- 
entations is about O.GlfcsT/a^, which corresponds to a 
TurnbuU coefflcicnt of 0.51. This average value is about 
10% higher than the value determined from nucleation 
rates on colloidal crystals, consistent with the differences 




FIG. 4. Step 4: Removing the cleaving walls. Pressure on 
the cleaving walls as a function of wall position for the three 
orientations of the interface. Filled and open symbols indicate 
pressure on the walls restraining fluid and crystal parts of the 
system, respectively. The error bars are of the order of the 
size of the symbols. 



TABLE I. Values of Wn for the steps of the thermodynamic 
integration process together with their sum 7. The statistical 
error estimates represent 2ct error bounds. 





(100) 


(110) 


(111) 




0.850 ± 0.001 


1.287 ±0.001 


1.125 ±0.001 




1.561 ±0.008 


1.989 ±0.007 


1.768 ± 0.008 


W4 


-1.789 ±0.005 


-2.639 ±0.006 


-2.311 ±0.005 


7 


0.62 ±0.01 


0.64 ±0.01 


0.58 ±0.01 



also about 10% higher than that determined from nucle- 
ation experiments on colloidal suspensions of polystyrene 
spheres, giving a rare comparison between direct evalu- 
ations of 7 and less accurate indirect determinations via 
nucleation theory. 
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found in other materials, such as bismuth (discussed 
above). Note that the hard-sphere 7 values are slightly 
anisotropic and increase in the order of (111), (100) and 
(110). That (111) has the lowest interfacial free energy 
is perhaps not surprising, since the (111) crystal face re- 
sembles most the structure adopted by the fluid against 
a structureless wall . 

It is generally accepted that the structure and ther- 
modynamics of dense simple fluids is dominated by the 
repulsive part of the potential, which can often be well 
approximated as a hard sphere, and the effect of the at- 
tractive part of the potential can be viewed as a small 
perturbation. If one considers the truncated Lennard- 
Jones system studied by Broughton and Gilmer and cal- 
culates an effective hard-sphere diameter at the triple 
point [T — 0.617e//cB), using the Barker-Henderson ap- 
proach j2^, one obtains a value of (0.39e/tT^) simply by 
rescaling the hard-sphere result calculated here. Thus, 
the attractive part of this potential accounts for only 
about 10% of the total 7, which is consistent with a sim- 
ilar observation by Curtin on the basis of DFT calcula- 
tions |l|l. 

To summarize, we have determined the crystal/melt 
interfacial free energy, 7 for the hard-sphere system di- 
rectly from simulation using a method that is similar to 
that Broughton and Gilmer ||l^ used for the truncated 
Lennard-Jones system except that we have replaced their 
external cleaving potentials with specially constructed 
cleaving walls. Although the method of cleaving walls 
is especially advantageous for the hard-sphere system, 
it could also be easily applied in modified form to con- 
tinuous potentials. The hard-sphere 7 obtained is only 
slightly dependent upon orientation and averages about 
O-GI/cbT/ct, consistent with some previous theoretical 
predictions from density-functional theory. This value is 
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